clear;clc;

load('Result_Baseline.mat');
%% Setup
% Information Structure
Info            =   struct();
Info.ParamList  =   {'RHO_M','Taylor_Pi','Taylor_ir'};
Info.ParamNum   =   length(Info.ParamList);
Info.StdVec     =   [1 -1/4]'*0.01;
Info.ValueList  =   struct('RHO_M',(0.05:0.1:0.95)', ...
                           'Taylor_Pi',(1.00:0.10:5.0)',...
                           'Taylor_ir',(0:0.1:0.9)');
Info.FixedParam =   struct('RHO_M',0.63,'Taylor_Pi',1.1,'Taylor_ir',0.90);
Info.Num        =   0;
UnitValue       =   zeros(1,Info.ParamNum);
for ii=1:Info.ParamNum
    vv              =   Info.ParamList{ii};
    Info.Idx.(vv)   =   Info.Num+(1:1:length(Info.ValueList.(vv)))';
    Info.Num        =   Info.Num+length(Info.ValueList.(vv));
    UnitValue(ii)   =   Info.FixedParam.(vv);
end

Info.BenchMark  =   UnitValue;
ParamMat        =   repmat(UnitValue,[Info.Num,1]);

for ii=1:Info.ParamNum
    vv              =   Info.ParamList{ii};
    ParamMat(Info.Idx.(vv),ii) ...
                    =   Info.ValueList.(vv);
end

%% Computation
TempFun         =   @(Param)SubFun_Param2Moment(PP,SS,MODEL,EquJac,Info,Param);

Result          =   cell(Info.Num,1);

parfor ii=1:Info.Num
    Result{ii}      =   TempFun(ParamMat(ii,:));
end

Result_Benchmark=   TempFun(Info.BenchMark);

save('Identification.mat','Info','Result','Result_Benchmark');

%% Collect the Results
ParamLabel      =   struct('RHO_M','$\rho_{M}$',...
                           'Taylor_Pi','$\lambda_{\pi}$',...
                           'Taylor_ir','$\lambda_{i}$');
StatList        =   {'C to $\varepsilon_{M}$','ir to $\varepsilon_{M}$',...
                     'CPI to $\varepsilon_{M}$'};
TempFun         =   @(IRF)[SubFun_MaxAbsVal(IRF.Eps_Eta_M.C)*100;...
                           SubFun_MaxAbsVal(IRF.Eps_Eta_M.ir)*400;...
                           SubFun_MaxAbsVal(cumsum(IRF.Eps_Eta_M.Pir))*100 ...
                           ];
StatMat         =   NaN(3,Info.Num);
for ii=1:Info.Num
    if ~isempty(Result{ii})
        StatMat(:,ii)   =   TempFun(Result{ii});
    end
    
end 
BenchMarkStat   =   TempFun(Result_Benchmark);
%% Figures
Folder          =   'TableGraphs/Identification/';
mkdir(Folder);
ParamList       =   fieldnames(ParamLabel);
N_col           =   length(ParamList);
N_row           =   size(StatMat,1);
fig_size        =   [N_col,N_row].*[1,1]*1/4; 
fig_gap         =   [0.03,0.03];
fig_VMargin     =   [0.05,0.05];
fig_HMargin     =   [0.05,0.05]; 

fig_SubPlot     =   @(ii_vv)subtightplot(N_row,N_col,ii_vv,...
                                         fig_gap,fig_VMargin,fig_HMargin);
Fig             =   FigureSetup(['Identification'],fig_size);
for ii=1:N_col
    vv          =   ParamList{ii};
    for jj=1:N_row
        ax          =   fig_SubPlot((jj-1)*N_col+ii);
        XAxis       =   Info.ValueList.(vv)';
        YAxis       =   StatMat(jj,Info.Idx.(vv));
        LineStyle   =   {'-';'b';2};
        plot(XAxis,YAxis,'-b','LineWidth',2);
        yline(BenchMarkStat(jj),':k','LineWidth',1);
        xline(Info.FixedParam.(vv),':k','LineWidth',1);
        ylabel(StatList{jj},'Interpreter','Latex','FontSize',8); 
        title(ParamLabel.(vv),'Interpreter','Latex','FontSize',8);
        xlim([min(XAxis),max(XAxis)]);
        ylim([min(min(StatMat(jj,:)),BenchMarkStat(jj)),max(max(StatMat(jj,:)),BenchMarkStat(jj))]);
        ax.FontSize =   8;
        ax.XGrid    =   'on'; ax.YGrid  =   'on';
%         ax.XAxis.Exponent   =   0;
%       ax.YAxis.Exponent   =   0;
%         ax.YAxis.TickLabelFormat    =   '%g2';
    end
end
print('-depsc','-r1000',Fig,[Folder,'Identification']);

